#############################################################################
# FE per period
# 3 periods:1950-1968, 1969-1987, 1988-2005
#
############################################################################

rm(list=ls())
graphics.off()

library(data.table)
library(ggplot2)
library(haven)
library(xtable)
library(fixest)
library(mvtnorm)
library(Hmisc)
library("RColorBrewer")

memory.limit(size=200000)

#############################################################################
# 1. Data and and FE where the slopes are constant within period
#############################################################################

# Data
load("Data_new/GFE/Data_for_groups.RData")
rm(dt_slop,dt_fe)
time0=Sys.time()

# Review basic FE regression
dim(dt_res)
dt_res[,lapply(.SD,wtd.mean,weight=corn_area)]
dt_res[,lapply(.SD,wtd.var,weight=corn_area)]

# Generate period variable
dt[year>=1950 & year<=1968, period:=1]
dt[year>=1969 & year<=1987, period:=2]
dt[year>=1988 & year<=2005, period:=3]
summary(dt$period)

# Generate interactions of period with ddr
dt[,ddr_p1:=ddr*(period==1)]
dt[,ddr_p2:=ddr*(period==2)]
dt[,ddr_p3:=ddr*(period==3)]
summary(dt[,.(ddr_p1,ddr_p2,ddr_p3)])

# Re run FE with slopes and county FE by period ()
t0=Sys.time()
fe = feols(yield ~  precr + ddr_p1:factor(fips) + ddr_p2:factor(fips) + ddr_p3:factor(fips)  | factor(fips)^factor(period) +
                    factor(year)^factor(state) ,data=dt,cluster="state")
summary(fe)
t1=Sys.time()
ft4 = t1-t0
ft4
rm(t0,t1,ft4)

# Source auxiliary functions
source("Code/0. CODE Auxiliary.R")

# Save the slopes in data table: not every slopes in every period
dt_slopes = slopes_dt_period(f=fe,np=3)
summary(dt_slopes)
dim(dt_slopes)
table(dt_slopes$period)

# Check number of coefficients: number of slopes + 1 precipitation
nobs(fe)
dim(fe$coeftable)
dim(dt_slopes)
summary(fe$coeftable)

# Marginal effects
dim(dt_res)
dt_res = merge(dt_res,dt[,.(fips,year,period)],by=c("fips","year"),all=TRUE)
dim(dt_res)
dt_res = merge(dt_res,dt_slopes[,.(fips,period,coef)],by=c("fips","period"),all=TRUE)
dim(dt_res)
summary(dt_res)
dt_res[,mg_fep:=coef]
summary(dt_res)

rm(f3,slopes_dt)
time1=Sys.time()
save.image("Data_new/FE_period_results.RData")

#############################################################################
# 2. Characterize results
#############################################################################

rm(list=ls())
load("Data_new/FE_period_results.RData")
save_res="Results"

# 113 observations without mg, corresponding to 112 counties, for whom one period slope is not identified
# 112 of those counties have 1 of their observations only in 1 of the periods
# 1 county has only 1 observation in period 1 and 1 observation in period 3
dt_res[,mg4:=mg_fep]
summary(dt_res)
summary(dt_res[is.na(mg4)])
length(unique(dt_res[is.na(mg4),fips]))
fips_now= unique(dt_res[is.na(mg4),fips])
table(dt_res[fips %in% fips_now,.(fips,period)])
dt_res[fips==48347,.(year,period,mg4)]

# Generate the mg3 restricted to same data
dt_res[,mg3r:=mg3]
dt_res[is.na(mg4),mg3r:=NA]
summary(dt_res)

#####################################
# Descriptive statistics (weighted)
#####################################

# Per period (only common subsample)
dt_aux=dt_res[is.na(mg3r*mg4)==FALSE,.(fips,year,period,mg3r,mg4,corn_area)]
summary(dt_aux)
dim(dt_aux)
dt_aux[period==1,mg3_p1:=mg3r]
dt_aux[period==1,mg4_p1:=mg4]
dt_aux[period==2,mg3_p2:=mg3r]
dt_aux[period==2,mg4_p2:=mg4]
dt_aux[period==3,mg3_p3:=mg3r]
dt_aux[period==3,mg4_p3:=mg4]
summary(dt_aux)
table(dt_aux$period)
length(unique(dt_aux[period==1,fips]))
length(unique(dt_aux[period==2,fips]))
length(unique(dt_aux[period==3,fips]))
nn = c("mg3_p1","mg3_p2","mg3_p3","mg4_p1","mg4_p2","mg4_p3")
s=as.matrix(dt_aux[,lapply(.SD,wtd.quantile,weights=corn_area,probs=c(0.1,0.25,0.5,0.75,0.9)),.SDcols=nn])
s=rbind(s,as.matrix(dt_aux[,lapply(.SD,wtd.mean,weights=corn_area),.SDcols=nn]))
s=rbind(s,as.matrix(dt_aux[,lapply(.SD,wtd.var,weights=corn_area),.SDcols=nn]))
rownames(s)=c(paste0("Percentile ",c(10,25,50,75,90)),"Mean","Variance")
rm(dt_aux)
s
dig=3
print(xtable(s,type="latex",digits=dig,display=c("s",rep("f",dim(s)[2]))),hline.after = NULL,
      file=paste0(save_res,"/Table_time_varying_FE.tex"),include.rownames = TRUE,include.colnames=FALSE,
      sanitize.text.function = function(x){x},only.contents = TRUE)

